Direct van der Waals simulation (DVS) of phase-transforming fluids

We present the method of direct van der Waals simulation (DVS) to study computationally flows with liquid-vapor phase transformations. Our approach is based on a discretization of the Navier-Stokes-Korteweg equations, which couple flow dynamics with van der Waals’ nonequilibrium thermodynamic theory of phase transformations, and opens an opportunity for first-principles simulation of a wide range of boiling and cavitating flows. The proposed algorithm enables unprecedented simulations of the Navier-Stokes-Korteweg equations involving cavitating flows at strongly under-critical conditions and 𝒪(105) Reynolds number. The proposed technique provides a pathway for a fundamental understanding of phase-transforming flows with multiple applications in science, engineering, and medicine.


INTRODUCTION
Flows of phase-transforming fluids are principal across science, engineering, and medicine. Management of electronics cooling, which depends heavily on liquid-vapor flows, remains a critical barrier to creating more powerful datacenter computers and meeting the performance demands of an increasingly computerized society and industry. The collapse of a cavitation bubble, which is another notable example of flows with phase transformations, has fascinated scientists for decades due to the extreme conditions generated, including temperatures of up to 5000 K, emission of light, and strong shock waves and jets (1,2). Although cavitation continues to be an important concern in the design of marine propellers, it has also been exploited technologically for ultrasonic cleaning and drug delivery (3,4). Despite their prevalence and importance, our understanding of fluid flows with phase transformations remains poor, partially due to the challenges that they pose to computational methods. Phasetransforming flows involve nonequilibrium thermodynamics, large viscosity, and density ratios, moving interfaces with topological changes and flow physics that spans a wide range of time and length scales. The most advanced computational methods are based on compressible flow models for mixtures of liquid and vapor. Although mixture models have been successful in several applications, their approach to phase change is either based on thermodynamic equilibrium or on phenomenological models that enter the mass balance equations and are known to have an important effect on the predictions (5). The latter phase change models, also called mass transfer functions, involve parameters that depend on the flow conditions and need frequent recalibration. These models cannot predict nucleation of vapor bubbles from pure liquid, which precludes further mechanistic understanding of, arguably, the most critical problem in cavitating and boiling flows (6).
Van der Waals proposed a first-principles thermodynamic theory of liquid-vapor phase change (7). The model is based on a nonconvex bulk Helmholtz free energy extended with a nonlocal term that accounts for interfacial energy. The use of a nonconvex bulk thermodynamic potential permits to incorporate the state-of-the-art theory of phase transformations that enables the prediction of nucleation and spinodal decomposition. Nonconvex potentials have found marked success in predicting thermodynamic properties and critical points of liquid-vapor mixtures (8). Van der Waals' thermodynamic theory can be coupled with the balance equations of compressible flows in a thermodynamically consistent manner that guarantees that the second law is satisfied for an arbitrary process compatible with the balance laws. The result of coupling van der Waals' theory with flow is the Navier-Stokes-Korteweg (NSK) equations. The potential of the NSK equations for mechanistic understanding and prediction of liquid-vapor flows has been exploited to study nucleation (9), fluid instability under shear (10), and bubble collapse (11). Although many numerical schemes have been proposed to solve the NSK equations, such as the local discontinuous Galerkin method (12,13) and relaxation models (14,15), current computational methods are limited to micrometer-scale flows without solid walls or flow conditions very close to criticality. Thus, the predictive capability of the NSK equations remains unrealized for a wide range of boiling and cavitating flows at length scales larger than a few micrometers.
Here, we present unprecedented three-dimensional simulations of wall-bounded cavitating flows at centimeter scale and O(10 5 ) Reynolds number using the NSK equations. Because our simulations are based only on van der Waals' thermodynamic theory and fundamental continuum mechanics without additional modeling assumptions, we call them direct van der Waals simulations (DVS). Our computations are enabled by a residual-based, stabilized discretization concept that does not require hyperbolicity of the isentropic form of the equations, extends to van der Waals fluids the Streamline Upwind Petrov-Galerkin (SUPG) technique (16,17) and the discontinuity-capturing (DC) operators (18), and improves the thickened interface methods (19,20). We illustrate the algorithm's performance with a parametric study of cavitating flow past a cylinder and a simulation of flow over a wedge that shows sheet-to-cloud transition. Our results are in good agreement with experiments, indicating that the proposed algorithm opens the opportunity to predict boiling and cavitating flows at centimeter scale or even larger using minimal modeling assumptions.

Model overview
The NSK equations are derived from the functional Helmholtz free energy Here, Ω is the fluid domain, ρ is the fluid's density, and ψ is the bulk Helmholtz free energy per unit volume, while λ and η are constants that control interfacial energy and interface thickness, respectively. The thermodynamic potential in Eq. 1 differs from standard potentials used for compressible flows in two critical aspects that are interconnected. First, ℋ depends not only on ρ but also on its gradient. Second, because ℋ depends on the density gradient, the system of equations remains well-posed even if ψ is nonconvex (21). The possibility of using a nonconvex bulk free energy per unit volume ψ allows us to use the state-of-the-art theory in nonequilibrium phase transformations. From the thermodynamic potential given in Eq. 1, we can derive the NSK equations using balance laws for mass, linear momentum, angular momentum, energy, and the second law of thermodynamics. The NSK equations for an isothermal system are ∂ρ ∂t þ r � ðρuÞ ¼ 0 ð2Þ where Eqs. 2 and 3 represent mass and linear momentum conservation, respectively. Here, u is the fluid velocity, p ¼ ρ 2 ∂ðψ=ρÞ ∂ρ is the fluid pressure, and I is the identity tensor. The tensor τ denotes viscous stresses, which, for a Newtonian fluid under Stokes' hypothesis, are given by where μðρÞ is the density-dependent viscosity coefficient; see Methods. The Korteweg stress tensor is and accounts for the interfacial stresses. The challenges in the simulation of Eqs. 2 to 5 for wall-bounded, large Reynolds number flows at centimeter scale emanate from two difficulties. First, there is a very large disparity between the length scale at which interfacial physics occurs and the largest length scale that controls flow physics. We address this by proposing the stabilized thickened interface method (sTIM); see Methods. Second, the inviscid NSK equations with vanishing Korteweg stress are not hyperbolic, which precludes the direct use of most standard computational methods for compressible flows. We bypass this difficulty using residual-based stabilization with shock capturing; see Methods.
To illustrate the potential of DVS, we study cavitating flow over a circular cylinder and over a wedge at centimeter scale. For all cases, we impose free-stream inlet boundary conditions (u ∞ and p ∞ ) using an acoustically absorbing sponge layer (22). The flow conditions are characterized by the free-stream cavitation number and Reynolds number.
The free-stream cavitation number is where ρ ∞ is the density that corresponds to p ∞ in our equation of state (EoS), and p v is the vapor pressure. The free-stream Reynolds number is R ' e ¼ ρ 1 u 1 '=μ l , where l is a problem-dependent length scale and μ l is the dynamic viscosity in the liquid phase. The strength and extent of cavitation will be measured using the void fraction α = (ρ ∞ − ρ)/(ρ ∞ − ρ v ).

Cavitating flow over a circular cylinder
A flowing fluid accelerates as it moves around the leading edge of a cylinder. The fluid's acceleration leads to a pressure drop that can trigger cavitation. Flows over cylinders have been often used to study cavitation because, depending on the free-stream conditions, they can feature different types of cavitation and different inception locations. Here, we perform a parametric study varying the freestream pressure to produce free-stream cavitation numbers that span the range σ ∞ = 0.25 (strong cavitation) to σ ∞ = 6.0 (no cavitation). Figure 1A shows snapshots of the instantaneous void fraction for different cavitation numbers under a flow field that goes from left to right. For σ ∞ = 3 (left panels), we observe cyclic cavitation. In this cavitation regime, the small cavities formed at the cylinder's surface detach almost instantaneously and are captured by the vortex immediately downstream of the cylinder. Because of their small sizes, these cavities collapse shortly after leaving the vortex. For transitional cavitation at σ ∞ = 1.25 (center column), some vapor pockets separate instantaneously from the cylinder's surface. Some cavities, however, remain attached to the cylinder for a time interval, grow, and eventually are carried downstream by the flow. Our simulation for σ ∞ = 0.25 shows fixed cavitation. In this case, a large fraction of the cylinder surface is consistently covered by vapor. The average cavity length remains stable over time, but its trailing edge continuously sheds gas pockets. The time-averaged vapor fraction 〈α〉 offers a more conclusive picture of the primary location of the cavity for each case; see Fig. 1B. For cyclic cavitation, the cavity is entirely detached from the cylinder. For transitional cavitation, the time-averaged cavity is attached to the cylinder and has a length that is comparable to the cylinder's diameter. For fixed cavitation, the cavity length is much larger than the cylinder and its thickness also exceeds the cylinder's diameter. Following (23), in Fig. 1C, we show the time-averaged length of the cavity 〈L〉 relative to the cylinder's diameter as a function of the ratio cavitation number σ * . The results are in good agreement with experiments (24) and past numerical studies (25). We observe that, although the cavity length decreases monotonically with the cavitation number in the majority of the plot, there is a small region close to the boundary between transitional and cyclic cavitation where it increases. This counterintuitive result has also been observed experimentally (24). On the basis of our results, one potential explanation is as follows: As the cavitation regime changes from cyclic to transitional, the size of the cavities attached to the cylinder grows. The presence of larger cavities at the cylinder's surface reduces the vortex strength and leads to weaker cavitation inside the vortex. Although for smaller cavitation number, larger cavities are shed into the free stream, they are short-lived because the free-stream pressure is relatively large and do not contribute notably to increase <L>. Thus, in this regime, the overall effect of the cavitation number increase is a larger cavity length. Figure 1D shows the point-wise, The computations are performed on a two-dimensional (2D) domain whose external boundary is an ellipse with a semimajor axis of 30D and a semiminor axis of 11.25D. The cylinder is located at the center of the ellipse. We use 77,274 C 1continuous quadratic elements to discretize the domain, where 243 elements are placed on the cylinder surface. An acoustically absorbing sponge layer with a width of 2D is placed near the edge of the ellipse. The temperature is T = 300 K. The free-stream velocity is u ∞ = 15.2 m/s, the cylinder's diameter is D = 2 mm, and the dynamic viscosity of the liquid phase is μ l ¼ 10 À 3 Pa·s, which corresponds to R D e ¼ 2:6 � 10 4 . The dynamic viscosity of the vapor phase is μ v ¼ 10 À 5 Pa·s; see Eq. 9. We vary the freestream pressure to change the cavitation number. We choose λ = 10 −16 m 7 /kg/s 2 and η = 10 7 as interfacial parameters, which yields the surface tension for a liquid-vapor interface in water at the problem's length scale. To reduce the computational cost, the initial condition is obtained from an incompressible flow simulation with density ρ ∞ , which implies that α = 0. (A) Instantaneous void fraction for free-stream cavitation number σ ∞ of 3.0 (cyclic), 1.25 (transitional), and 0.25 (fixed). (B) Time-averaged void fraction 〈α〉. (C) Average vapor cavity length nondimensionalized by cylinder diameter as a function of the ratio cavitation number σ * = (σ ∞ − σ ck )/(σ i − σ ck ), where σ ck and σ i are the choke and incipient cavitation numbers, respectively (23). Experimental data are obtained from (24). (D) Time-averaged local cavitation number 〈σ θ 〉 distribution on the cylinder. time-averaged cavitation number on the cylinder surface, Here, θ is a parametric coordinate along the cylinder's surface such that θ = 0 ∘ and θ = 180 ∘ correspond to the leading and trailing edges, respectively. For cyclic cavitation, 〈σ θ 〉 reaches a local minimum at θ ≈ 80 ∘ . For slightly larger values of θ, the pressure first increases due to flow deceleration and later decreases due to cavity shedding. For transitional cavitation, 〈σ θ 〉 decreases monotonically with θ, which reinforces the idea that cavitation inception is caused by instantaneous pressure fluctuations. In contrast with the previous two cases, for fixed cavitation, 〈σ θ 〉 drops abruptly to zero at θ ≈ 80 ∘ and remains at this value on the rest of the cylinder's surface. These results further emphasize the difference between the three cavitation modes. We observe that 〈σ θ 〉 has a sharp increase at θ ≈ 55 ∘ for fixed cavitation. To better understand this phenomenon, we show the time-averaged velocity magnitude for the entire cylinder (top) and near the separation point (bottom) in Fig. 2. The velocity inside the vapor pocket remains close to zero, which indicates that only a small fraction of the momentum is transported across the liquid-vapor interface. When we have a cavity consistently attached to the cylinder, more kinetic energy accumulates upstream and is converted into internal energy. Such conversion causes a local increment in pressure and a stronger adverse pressure gradient, which leads to the thickening of the boundary layer and earlier flow separation. Such a phenomenon has been observed experimentally (26,27) but has remained elusive for computational methods.

Sheet-to-cloud transition in cavitating flow over a wedge
Flows over a wedge have been often used to study cavitation problems. Figure 3A shows a schematic configuration of this physical system and our simulation setup. Under these conditions, the inlet flow accelerates along the wedge, which leads to a pressure drop that triggers cavitation. The cavity initially grows attached to the bottom wall, developing the shape of an elongated sheet. The sheet grows longer until it pinches off and transitions to a cloud. The cloud is a three-dimensional structure with features that range across multiple length scales. As the cloud travels downstream, it encounters increasingly large pressures that lead to bubble collapse, which generates jets and sound. The results are in agreement with the experimental observations (28) but reveal important aspects of the cavitation inception process and the sheet-tocloud transition. Understanding the flow conditions that trigger cavitation remains an outstanding challenge. Our results point to a complex scenario in which cavitation is a strongly unsteady and heterogeneous process that is tightly controlled by localized and instantaneous reductions of pressure. Figure 3B shows that the timeaveraged pressure remains well above the vapor pressure, but it is instantaneous descents of the pressure, at a level similar to the vapor pressure, that trigger cavitation. Figure 3B also illustrates that the instantaneous pressure decreases quickly along the wedge due to flow acceleration, but it does not reach a minimum at the wedge apex. Instead, the boundary layer separation that occurs downstream of the apex generates vortices that undergo stretching and further reduce the pressure, eventually leading to the formation of a vapor cavity. The transition from sheet to cloud cavitation is important because cavitation clouds have a higher potential to generate shock waves and noise. However, the mechanisms that control the transition remain poorly understood (29). Recent research (30,31) points to a scenario in which, as the cavity grows, the sheet becomes unstable and transitions into a cloud due to a combination of a reentrant jet and a condensation shock that travels upstream. Callenaere et al. (32) identify two transition types based on the sheet's thickness. In thick sheets, the jet plays a minor role until it reaches the cavity's leading edge and triggers the transition. In contrast, thin cavities break into smaller-scale, three-dimensional structures immediately after they are impinged by the reentrant jet. Figure 3C shows snapshots of the spanwise-averaged instantaneous void fraction. As the sheet cavity travels downstream, the reentrant jet starts to develop due to the presence of an adverse pressure gradient. Because the sheet is thick, it remains intact as the jet travels through. Once the jet fully penetrates the sheet, a cloud cavity pinches off the rest of the sheet. As the cloud cavity travels downstream, the remaining sheet starts to interact with free-stream nuclei and forms a second cloud cavity. Meanwhile, a thinner sheet cavity develops near the wedge apex due to pressure fluctuations. While this secondary sheet cavity develops, a reentrant jet is formed. Because the secondary sheet cavity is thinner, the reentrant jet immediately destabilizes it, leading to many smaller-scale three-dimensional structures.

DISCUSSION
We propose an algorithm that allows DVS of phase-transforming fluids for wall-bounded flows far from criticality and large Reynolds numbers at unprecedented length scales. Our algorithm is based on a residual-based formulation and a stabilized thickened interface method (sTIM). The proposed approach successfully addresses two critical challenges that limited existing computational methods, namely, the nonhyperbolic eigenstructure of the inviscid equations without Korteweg stress and the disparity of length scales between interfacial physics and flow physics. The strength of DVS is that it couples flow dynamics with a fundamental nonequilibrium theory of phase transformations without resorting to phenomenological approaches that require flow-dependent parameter calibration. DVS opens the possibility to gain mechanistic understanding of the most critical processes of phase-transforming flows, including nucleation of the vapor phase in boiling and cavitation.
To illustrate our approach, we performed a parametric study of flow over a circular cylinder, varying the free-stream pressure. As the free-stream pressure is reduced, DVS predicts a transition from noncavitating to cavitating flow. DVS also predicts the progression from cyclic to fixed cavitation in quantitative agreement with experiments. Our DVS results indicate that, as the vapor cavity attached to the cylinder's trailing edge grows larger, the separation point moves upstream. This subtle, yet critical phenomenon, has been observed in experiments but not in state-of-the-art cavitation simulations.
We performed a three-dimensional simulation of cavitating flow over a wedge of 1.5 cm height. Our DVS results capture a highly turbulent flow and the transition from sheet to cloud cavitation. The simulation shows that cavitation inception is tightly controlled by local pressure fluctuations and vortex dynamics. In agreement with experiments, DVS shows that thin and thick sheet cavities respond differently to reentrant jets and condensation shocks, which leads to distinctive destabilization mechanisms of the sheet cavity. Overall, our results highlight the predictive capabilities of Fig. 3. Sheet-to-cloud transition in cavitation over a wedge. The wedge height is H = 1.5 cm, and the flow temperature is T = 300 K. The dynamic viscosity of the liquid and vapor phases is μ l ¼ 10 À 3 Pa·s and μ v ¼ 10 À 5 Pa·s, respectively; see Eq. 9. The free-stream conditions are u ∞ = 15.17 m/s and p ∞ = 101,325 Pa, which correspond to R H e ¼ 2 � 10 5 and σ ∞ = 1. We choose λ = 10 −16 m 7 /kg per s 2 and η = 10 9 as interfacial parameters, which yields the surface tension for a liquid-vapor interface in water at the problem's length scale. To reduce the computational cost, the initial condition is obtained from an incompressible flow simulation with density ρ ∞ , which implies that DVS, which are particularly noteworthy because the modeling assumptions are minimal. We believe that DVS opens possibilities not only to simulate and predict flows of phase-transforming fluids but also to fundamentally understand bubble nucleation and cavitation inception.

Governing equations
The isothermal NSK equations can be written as Here, an inferior comma denotes partial differentiation (e.g., U ,t = ∂U/∂t) and repeated indices indicate summation over the spatial dimensions (e.g., F adv where x i denotes the ith Cartesian coordinate and d is the number of spatial dimensions). The vector U = [ρ, ρu 1 , ρu 2 , ρu 3 ] contains the conservation variables. The vectors F adv i , F diff i , and F c represent the advective fluxes, the diffusive fluxes, and the Korteweg stress, respectively, and they are defined as where δ ij is the Kronecker delta. In Eq. 8, we have used the identity ∇ · ζ = ληρ ∇ (Δρ). In the viscous stress tensor τ ¼ μðρÞ ru þ r T u À 2 3 r � u I À � , the dynamic viscosity is defined as where μ l=v and ρ l/v are the dynamic viscosity and saturation density for the liquid and vapor phases, respectively. To derive our algorithm, we define the primitive variables Y = [ρ, u 1 , u 2 , u 3 ] and the following transformation matrices whose explicit expressions are given in the Supplementary Materials. Using the transformation matrices, we can rewrite Eq. 6 in quasi-linear form

Cubic EoS
Cubic EoS are widely used to represent liquid-vapor equilibrium (8). The first cubic EoS is due to van der Waals (7), but many variants and extensions have been proposed thereafter, including the Soave-Redlich-Kwong (33) and Peng-Robinson models (34).
Here, we use the EoS where R is the specific gas constant, T is the temperature that is a constant for isothermal conditions, and a ( where T R = T/T c , the critical temperature is T c = 647.1 K, and

Stabilized thickened interface method
For temperatures below the critical temperature, equilibrium solutions of the NSK equations with Eq. 15 predict a liquid-vapor interface described by a continuous variation of density. At room temperature, the model predicts an interface thickness of less than 100 nm, in agreement with experiments and molecular dynamics simulations (36,37). In a simulation of the NSK equations, the interface thickness needs to be resolved by the computational mesh that implies that a three-dimensional centimeter-scale computation would require at least ∼10 16 degrees of freedom that is prohibitive in today's computer architectures. Enlargement of the interface can be achieved by increasing the parameter η in the governing equations. However, increasing η without modifying the EoS leads to an overprediction of surface tension that would make the results invalid.
Notably, the use of the thickened interface method (19,20) permits to enlarge the interface thickness while keeping surface tension constant. This is accomplished by increasing η and modifying accordingly the binodal region of the EoS. While the thickened interface method opens the possibility to perform larger-scale computations, it leads to the use of a nondifferentiable EoS. The lack of smoothness in the EoS leads to the formation of strong spurious shock waves at the interface that propagate throughout the computational domain and become a source of instability. To address this issue, we propose the sTIM. The formulation of sTIM is where A v=l ¼ ξ 1À η η ∂p EoS ∂ρ ðρ v=l ; TÞ, and ρ v and ρ l are the vapor and liquid saturation densities at temperature T, respectively. Equation 18 shows that, in the binodal region, ρ v < ρ < ρ l , the pressure is modified using the approach proposed in (20) and, thus, by increasing η, one can enlarge the interface while keeping surface tension constant. In the vapor (0 < ρ ≤ ρ v ) and liquid (ρ ≥ ρ l ) phases, the original EoS is modified with a stabilizing term whose strength is controlled by the parameter ξ. The stabilizing term is designed such that the following conditions are satisfied Equations 19 and 20 guarantee that saturation pressure remains unchanged, while Eqs. 21 and 22 ensure that the pressure is a differentiable function at saturation conditions for all ξ ≠ 0. When the stabilizing term is absent (ξ = 0), the sTIM reduces to the methodology proposed in (20). In our computations, we took ξ = 0.01, which guarantees that the EoS is smooth and produces changes in the pressure outside of the binodal region that are negligible. Figure 4 shows a plot of p and p EoS as functions of the density for ξ = 0.01 and several values of η. We can see that outside of the binodal region p is indistinguishable from p EoS . In the binodal region, p is different from p EoS for η ≠ 1 to achieve the desired effect of decoupled interface thickness and surface tension. The larger is η, the flatter is p in the binodal region. Our proposed sTIM reconstruction has two additional important properties: (i) The saturation conditions of the original EoS are exactly preserved and (ii) the reconstructed EoS still satisfies the Maxwell's rule of equal area exactly.

Variational operators Galerkin operator
The proposed computational method is based on a weak form of the NSK equations that is stabilized with residual-based terms. Our weak formulation makes use of several semilinear forms. The first one, which emanates from the weak form of Eq. 14 without stabilizing terms, is defined as where W ∈ V is a vector-valued weight function, V is a suitably chosen functional space, Ω is the computational domain, Γ is the boundary of Ω, and n i is the ith Cartesian coordinate of the unit outward normal to Γ.

SUPG operator
Stremline-Upwind/Petrov-Galerkin (SUPG) is a finite element stabilization method for advection-dominated flow that is applicable to incompressible and compressible flows (38,39). SUPG is a residual-based stabilizing scheme that provides stable solutions retaining optimal rate of convergence. The analysis performed in (40) provides a strong mathematical foundation for SUPG. Let us assume that the domain Ω is divided into N el elements each denoted by Ω e . We define the SUPG operator as Here, is the residual of the governing equations, τ SUPG ¼ A À 1 0τ SUPG is the stabilizing matrix for the primitive variables, andτ SUPG is the stabilizing matrix for the conservation variables, which is defined as (16,17) In Eq. 26, Δt is the time step size, C I is a positive constant derived from an element-wise inverse estimate (41), and G ij represents the components of the element metric tensor G, that is where x(Ξ) is the element isoparametric mapping. The matricesÂ � i andK ij are the conservation variable counterpart of A � i and K ij , which can be obtained aŝ Because the SUPG operator defined in Eq. 24 is residual-based, the matrices A � i can be chosen in multiple ways without compromising the accuracy of the algorithm. However, a poor choice of A � i will have a detrimental effect on the stability of the scheme. In classical gas dynamics, this choice is guided by an eigenvalue analysis of the isentropic system. The isentropic NSK equations are not hyperbolic on the entire phase space and the eigenvalue analysis cannot be used. On the basis of scaling arguments and local equilibrium at the liquid-vapor interface, we choose where A c;eq i is an approximation to A c i . In Eqs. 23 to 25, A c i takes on the form where e i is the ith vector of the Cartesian basis in dimension d + 1. The computation of A c i is ill-conditioned, especially in the bulk phases because both numerator and denominator approach zero. Using A c i to compute the A � i matrices leads to small perturbations in the numerical solution that are eventually amplified unless the time step is extremely small. To derive an approximation to A c i , we proceed as follows: Under equilibrium conditions, the equation p ,i − ληρΔρ ,i = 0 is satisfied. Doing basic manipulations, one can show that p ,ρ = ληρΔρ ,i /ρ ,i , where no sum on i is implied. In addition, we know that, under equilibrium, p ,ρ ≤ 0 in the interfacial region and p ,ρ ≈ 0 in the bulk phase. Thus, we define where p À ;ρ ¼ minð0; p ;ρ Þ. Unlike A c i in Eq. 30, the matrix A c;eq i in Eq. 31 can be computed in a stable manner at all points in the computational domain. Although our derivation of the matrices A � i assumes that the interface is under local equilibrium conditions, this assumption does not compromise the accuracy of the algorithm in any way. Our discretization method still features high-order accuracy because the A � i matrices are used only in the SUPG operator that also involves the residual.
This completes the definition of all the matrices on the righthand side of Eq. 26. To calculateτ SUPG , we need to compute the square root of a (d + 1) × (d + 1) matrix. In our simulations, this is done using the Denman-Beavers algorithm (42,43). Discontinuity capturing While the use of SUPG ensures stability and accuracy when the solution is smooth, it does not resolve effectively flow fields with sharp layers (44). To stabilize solutions with sharp layers, the SUPG formulation is usually augmented by a residual-based discontinuity capturing (DC) operator. When the solution is under-resolved by the mesh, the term adds local dissipation along direction of the solution gradient (45), which enhances the stability of the numerical method while retaining optimal rate of convergence (46). The DC  operator for primitive variables is given by (17,47) Here,κ DC ¼κ C e i � e i þκ M e iþ1 � e iþ1 is a (d + 1) × (d + 1) diagonal matrix with entrieŝ and C C and C M are O(1) positive constants for which we used the value C C = C M = 0.1. In Eqs. 35 and 36, |Res 1 (Y)| is the absolute value of the residual of the mass conservation equations, ‖ · ‖ denotes the Euclidean norm of a vector, Res 2:d+1 (Y) is the residual of the linear momentum balance equation, p þ ;ρ ¼ max ðp ;ρ ; 0Þ, and u rel = u − u ∞ is the relative velocity with respect to the free-stream velocity. DC schemes were developed in the context of classical gas dynamics. In liquid-vapor flows, however, the changes in density and viscosity are much more marked, and this has to be considered while designing the DC operator. For this reason, we proposed a scaling term β, which is designed to minimize numerical dissipation in the liquid phase while retaining stability in the vapor phase. We use the expression where ρ m represents a very small value of the density and β max is a constant that sets the maximum strength of the DC operator. In our simulations, we take ρ m = 0.01 kg/m 3 and β max = 1000. In the liquid phase, the speed of sound is high, the solution is primarily smooth, and the use of the DC operator is not necessary. In the interfacial region (ρ v < ρ < ρ l ), the value of β varies linearly between zero and one. In the vapor phase (ρ ≤ ρ v ), the fluid is highly compressible and the use of a robust DC is necessary to retain numerical stability. For densities in the range ρ m < ρ ≤ ρ v , we set β = 1, which is a commonly used value in gas dynamics simulations (18). The DC is maximum when ρ ≤ ρ m to avoid the appearance of negative densities.

Fully discrete formulation
The NSK equations include third-order derivatives of the density. Thus, for the operators introduced in Eq. 23 to be well defined, we need a discrete functional space that is at least globally C 1 -continuous. Although we use a spatial discretization based on isogeometric analysis that offers this capability (48), classical finite elements do not support globally C 1 -continuous spaces on complex three-dimensional geometries. Thus, to make our algorithm applicable to classical finite elements, we use the split approach that is based on introducing the additional unknown By treating μ as an independent unknown, we can redefine the matrix A c i as and rewrite the NSK equations as a larger system of equations with derivatives of order less or equal than two. To formulate our semidiscrete problem, we use a finite element space V h that satisfies the Dirichlet boundary conditions and V h 0 , an analogous discrete space that satisfies homogeneous conditions at the Dirichlet boundary. The semi-discretized problem is as follows: We used the generalized-α method to perform time integration (49). At each time step, the nonlinear system of equations is solved using the Newton-Raphson's method with a relative tolerance of 2.5 × 10 −4 . The linear systems of equations are solved using the generalized minimal residual method (50) with an additive Schwarz preconditioner. The time step is varied throughout the simulation to achieve convergence of the Newton-Raphson algorithm in three to four iterations. Our code makes use of the open-source package PETSc (51) and PetIGA (52).

Stability and accuracy of the proposed algorithm
Because the SUPG operator in Eq. 24 vanishes when the residual Res(Y) is zero, the matrices A � i in Eq. 24 can be chosen in multiple ways without compromising the rate of convergence of the algorithm. However, a poor choice of A � i can make the algorithm unstable for mesh sizes or time steps that are not sufficiently small to reach the asymptotic regime of the algorithm. Here, we show that two choices of A � i that are logical extensions of the matrices used in standard gas dynamics simulations render unsatisfactory results, while our choice, given by Eq. 29, produces vastly superior results. The alternatives to Eq. 29 that we study here are the following: (a) Case (a) corresponds to the standard SUPG operator used for compressible Navier-Stokes; see (17). Case (b) represents a plausible, but unsuccessful, extension of the SUPG operators from compressible Navier-Stokes to the NSK equations.
To compare these three algorithms, we simulate the dynamics of three vapor bubbles; see Fig. 5. The expected solution is that the bubbles will collapse one after another from the smallest to the largest. As shown in Fig. 5B (left column), when we use method (a), the two smallest bubbles collapse, but the largest bubble acquires an irregular shape and periodically oscillates until the simulation becomes unstable. Algorithm (b) produces results that, at the scale of the plot, are indistinguishable from those of the proposed algorithm (see center and right columns of Fig. 5B). However, the time step required to get convergence of the Newton-Raphson scheme is ∼200 times smaller when we use algorithm (b) than when we use the proposed method; see Fig. 5C. Overall, this shows that the proposed SUPG operator vastly outperforms naive extensions of the classical SUPG method to the NSK equations.
We now perform an additional numerical test to evaluate the numerical dissipation introduced by our algorithm. Most successful algorithms for compressible flows introduce numerical dissipation. However, to obtain an accurate method, the amount of numerical dissipation should quickly approach zero as the mesh is refined. We study the artificial dissipation of our method by simulating the oscillation of an inviscid, planar liquid-vapor interface driven by an initial disturbance. Because the flow is inviscid, we expect a periodic oscillation of the interface without any decay in the amplitude of the disturbance. The initial velocity is zero and the initial void fraction is depicted in Fig. 6A. We show snapshots of the void fraction in the region of interest at multiple times in Fig. 6B. These pictures show a periodic oscillation of the interface. A more informative description of the oscillation is given by Fig. 6C (left), which shows the time evolution of α at the center of the domain for different mesh sizes. The plot shows that the simulation remains stable even for extremely coarse meshes (N = 32), but there is a notable decay in the wave amplitude due to numerical dissipation. The amplitude decay is less noticeable as we refine the mesh and becomes very small for our finest mesh (N = 256), even after 10 complete wave periods.